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We investigate the interplay between frustration and zero-point quantum fluctuations in 
the ground state of the triangular and J1 — J2 Heisenberg antiferromagnets, using finite- 
size spin-wave theory, exact diagonalization, and quantum Monte Carlo methods. In the 
triangular Heisenberg antiferromagnet, by performing a systematic size-scaling analysis, 
we have obtained strong evidences for a gapless spectrum and a finite value of the 
thermodynamic order parameter, thus confirming the existence of long-range Neel order. 
The good agreement between the finite-size spin-wave results and the exact and quantum 
Monte Carlo data also supports the reliability of the spin-wave expansion to describe 
both the ground state and the low-energy spin excitations of the triangular Heisenberg 
antiferromagnet. In the J1—J2 Heisenberg model, our results indicate the opening of 
a finite gap in the thermodynamic excitation spectrum at J2/J1 — 0.4, marking the 
melting of the antiferromagnetic Neel order and the onset of a non-magnetic ground 
state. In order to characterize the nature of the latter quantum-disordered phase we 
have computed the susceptibilities for the most important crystal symmetry breaking 
operators. In the ordered phase the effectiveness of the spin-wave theory in reproducing 
the low-energy excitation spectrum suggests that the uniform spin susceptibility of the 
model is very close to the linear spin-wave prediction. 

Keywords: The contents of the keywords 

1. Introduction 

The physics of quantum antiferromagnets is a very old topic, dating back to the 
early days of quantum mechanics itself. Nonetheless, after many years of intensive 
study, the interest in this research field is still high, with several new problems aris- 
ing from the behavior of low-dimensional magnetic materials. This is also due to the 
existence of simple toy-models in which the interplay between antiferromagnetism, 
symmetry, dimensionality and strong quantum correlations leads to fascinating ef- 
fects in the low-temperature physics, often reproducing the behavior of real systems. 
Among them, the nearest-neighbor Heisenberg Hamiltonian, 




(1) 



n.n. 



where 



!f,S?,S?) are spin-s operators and J is the (positive) exchange inte- 
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gral, has certainly played a central role as an ideal test ground to investigate the 
influence of quantum effects on the mechanism of spontaneous symmetry break- 
ing. In fact, in contrast to the ferromagnetic state, the classical Neel state is not 
an eigenstate of the Heisenberg Hamiltonian and therefore, in general, the ground 
state of the latter does not have a purely classical representation. Hence, quantum 
effects may play an important role in modifying the zero-temperature properties of 
the model from the classical (s — > oo) limit. In particular, reduced dimensionality 
and a small spin value might enhance zero-point quantum fluctuations up to the 
point of destroying the classical Neel order, thus stabilizing a ground state with 
symmetries and correlations different from its classical counterpart. 

Indeed, in one dimension and for s = 1/2, a famous exact solution found by 
Bethe in 1931 (Ref. Ill) showed that quantum effects prevent the onset of true long- 
range antiferromagnetic order, giving instead a power-law decay of the spin-spin 
correlation functions. Despite Bethe's promise to generalize his solution to the 
two-dimensional square lattice case, appearing in the conclusions of his paper, this 
was never done, and the issue of the existence of long-range order in the ground 
state of the two-dimensional Heisenberg model has been left unsolved for many 
years. The rigorous proof of the ordered nature of the ground state of the square 
Heisenberg antiferromagnet was given in fact, for s > 1, only in 1986 J3 and has not 
been extended yet to the spin-half case where zero-point quantum fluctuations are 
stronger. 

This problem became a hot topic when possible connections between a non- 
magnetic ground state and the mechanism of high-T c superconductivity were put 
forward by Anderson in 1987.0 In fact, since the stoichiometric compounds of the 
high-T c superconductors are good realizations of a s — 1/2 square Heisenberg anti- 
ferromagnet, this conjecture focused the attention on the properties of this system. 
Fortunately enough, at that time the development of modern computers was such 
that the use of numerical techniques could compensate for the lack of exact ana- 
lytical results. In particular, quantum Monte Carlo methods have been of crucial 
importance, by allowing one to perform a systematic size-scaling of the physical 
observables and therefore to reach a definite conclusions As a result, even if a rig- 
orous proof is still lacking, there is at present a general consensus about the ordered 
nature of the ground state of the spin-half square Heisenberg antiferromagnet: in 
two dimensions, reduced dimensionality and a low spin value do not seem enough 
to stabilize, within the Heisenberg model, a non-magnetic ground state. 

Better candidates for a realization of disordered ground states in two dimensions 
are frustrated spin models. In these systems, in fact, the usual antiferromagnetic 
alignment between spins is hindered by the geometry of the lattice or by the presence 
of competing interactions. As a result, a general feature introduced by frustration 
is a less stable classical minimum energy configuration which is more likely to be 
destabilized by zero-point quantum fluctuations for a small spin value. Among this 
class of systems two prototypical examples are given by the triangular Heisenberg 
antiferromagnet, and the J\ — Ji Heisenberg model. The nature of the ground state 
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Figure 1: The classical Neel state consists of coplanar spins forming ±27r/3 angles 
between nearest neighbors. This leads to a V3 x V3 periodicity with the spins on 
the three sublattices A,B,C ferromagnetically aligned. 



in these frustrated spin models represents the main topic of this paper. 

The triangular Heisenberg antiferromagnet is described by the Hamiltonian (Q) , 
where i and j are the sites of a triangular lattice. Due to the geometry of the lattice 
(see Fig. Q), the classical minimum energy configuration of this model is not the 
usual Neel state with antiparallel spins on neighboring sites. In fact, if two spins 
on an elementary triangular plaquette minimize their exchange energy by aligning 
antiparallel, the third one cannot do the same because it cannot be antiparallel 
to both of them, simultaneously. As a result, the minimum energy configuration 
consists of coplanar spins forming ±27r/3 angles between nearest- neighbors and this 
leads to a x -\/3 periodic Neel state with the spins ferromagnetically aligned on 
each of the three sublattices (Fig. |l|). The resulting state, having an energy per 
bond twice than the optimal one, is far less stable than that on the square lattice. 

In the J1—J2 model, instead, frustration arises on the square lattice because of 
the presence of competing interactions, the Hamiltonian being 

H = Ji J2 ^ ■ Sj + J 2 J2 S ^ , ( 2 ) 

n.n. n.n.n. 

where J\ and J2 are the antiferromagnetic couplings between nearest- and next- 
nearest- neighbors, respectively. Classically, the minimum energy configuration has 
the conventional Neel order for J2/J1 < 0.5 [Fig. || (a)]. By increasing further the 
frustrating interaction J2 this configuration is destabilized and, for J2/J1 > 0.5, the 
system decouples into two Neel ordered sublattices. At the purely classical level, 
the energy of the latter configuration is independent of the relative orientations 
of the staggered magnetizations on the two sublattices. However, this degeneracy 
is partially lifted by zero-point quantum fluctuations even at the lowest order in 
1/s so that in the s — > 00 limit the minimum energy configuration is the so-called 
collinear state [Fig. |^ (b)] with the spin ferromagnetically aligned in one direction 
and antiferromagnetically in the other, corresponding to a magnetic wavevector 
Q = (tt, 0) or Q = (0, 7r)Jj Exactly at J2/J1 = 0.5 any classical state having zero 
total spin on each elementary square plaquette is a minimum of the total energy. 
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Figure 2: The two sublattices Neel (a) and the collinear (b) classical states. 



These states include both the Neel and the collinear states but also many others 
with no long-range order. The occurrence of a non-magnetic ground state in the 
quantum case, for a small spin value, is therefore likely around this value of the 
J2/J1 ratio. 

The recent experimental finding of real compounds described by the triangular 
and the J\—Ji Heisenberg antiferromagnets have renewed the interest in these frus- 
trated spin systems. In particular, the K/Si(lll)-\/3 x v^-B interfaced has turned 
out to be a good experimental realization of a spin-half Heisenberg antiferromagnet 
on a triangular lattice. In fact, due to strong electronic correlations, the surfaces 
states consist of a triangular arrangement of half-filled dangling bonds, which are 
localized and carry local s = 1/2 magnetic moments coupled antiferromagnetically. 
Recent experimental realizations of the spin-half J\ — Ji Heisenberg model have been 
found instead in the Li2VOSi04 and Li2VOGe04 compounds.El These are three- 
dimensional systems formed by stacked square planes of V (s — 1/2) ions with a 
weak inter-plane interaction. The structure of the V planes suggests that both 
the superexchange couplings between first and second neighbors can be significant 
and indeed, the first experimental results have indicated that these antiferromag- 
netic couplings are of the same order of magnitude.El In addition, the possibility of 
performing measurements under pressure will also allow one in the near futureO to 
tune the J2/J1 ratio and to investigate the properties of these systems in various 
regimes of frustration. 

In this work the problem of the nature of the ground state of these frustrated 
spin systems is tackled using various techniques, namely: the finite-size spin-wave 
theory, exact diagonalization of small clusters by the Lanczos algorithm, and several 
zero-temperature quantum Monte Carlo methods. The finite-size spin-wave theory 
has been recently proposedli!] as a generalization of one of the oldest analytical 
techniques in the study of quantum magnetism.til In particular, this spin-wave ex- 
pansion allows one to deal with finite clusters while avoiding the spurious Goldstone 
modes divergences in a straightforward way. Even if these results are biased by the 
long-range order hypothesis, nevertheless useful information on the thermodynamic 
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ground state can be extracted from the comparison with the numerical results on 
finite systems or from the occurrence of a breakdown of the 1/s expansion. 

The Lanczos methodEl allows the exact evaluation of static and dynamical prop- 
erties of the finite-size system and, especially when combined with a careful analysis 
of the symmetry of the low-energy excited statesjlj can provide clear indications 
about the nature of the ground state. However, due to memory constraints, ex- 
act diagonalizations techniques are limited in two dimensions to very small clusters 
(~ 30 sites) so that it is in general difficult to perform a systematic size scaling 
of the important physical observables. In order to numerically investigate larger 
systems, different approaches are therefore necessary. 

In the unfrustrated cases, quantum Monte Carlo has turned out to be an essential 
instrument for studying both thp_ground-state and the finitc-tcmperature properties 
of a quantum antifcrromagnet.E-H Unfortunately, in the frustrated cases standard 
stochastic techniques cannot be applied, as their reliability is strongly limited by the 
well-known sign problem. This numerical instability originates from the vanishing of 
the signal-to-noise ratio in the Monte Carlo sampling which occurs within bosonic 
models in the presence of frustration or, in general, in fermionic systems. 

Presently, the sign problem can be controlled only at the price of introduc- 
ing some kind of approximation. Apart from purely variational calculations, the 
simplest approximation scheme in the framework of one of the most efficient zero- 
temperature algorithms - the Green function Monte CarlcEj - is the fixed-node 
(FN) technique.Ej In this technique, the exact imaginary time propagator e~ Tn - 
used to filter out the ground state from the best variational guess \ipc) ~ is replaced 
by an approximate propagator e~ T ^ FN such that the nodes of the propagated state 
e T ^ FN IV^g) do not change, due to an appropriate choice of the effective FN Hamil- 
tonian 7Yfn (which in turn depends on IV'g))- The FN approximation becomes 
exact if the so called guiding wavefunction IV'g) is the exact ground state. However, 
the fixed-node results are usually strongly biased by this ansatz so that it is in 
general difficult to extract reliable information about the ground-state correlations 
whenever they are not well reproduced by the variational guess. 

In order to overcome this difficulty, here we have used the recently developed 
Green function Monte Carlo with Stochastic Reconfiguration (GFMCSR),lHI'tL3 which 
allows one to release the fixed-node approximation in a controlled way and to obtain 
much more accurate estimates of the ground-state correlations, thus reproducing 
also ground-state properties that are not contained at the variational level in the 
guiding wavefunction. During each short imaginary time evolution r — > r + At, 
where both the exact and the approximate propagation can be performed with- 
out sign problem instabilities, the FN dynamic i s sys tematically improved by re- 
quiring that a given number p of mixed averagestJx3 of correlation functions are 
propagated consistently with the exact dynamic. By increasing the number of cor- 
relation functions one typically improves the accuracy of the calculation since the 
method becomes exact if all the independent correla tion functions are included in 
the stochastic reconfiguration (SR) scheme. Typically,E3El few correlation functions 
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(p ~ 10) allow to obtain rather accurate values of the ground-state energy with an 
error much less than 1% - for lattice sizes (N < 36) where the exact solution is 
available numerically - and without a sizable loss of accuracy with increasing size. 
Such accuracy is usually enough to reproduce some physical features that are not 
contained at the variational level. 

In this paper, using these numerical techniques, after a brief review of the basic 
concepts of spontaneous symmetry breaking in a quantum antiferromagnet (Sec. 2), 
we provide two clear examples of systems in which the combined effect of frustra- 
tion and quantum fluctuations do not or do change the zero-temperature long-range 
properties of their classical counterparts. In particular, we find that the thermody- 
namic ground state of the spin-half triangular Heisenberg antiferromagnet (Sec. 3) 
is most likely long-range ordered although with a remarkable reduction of the order 
parameter with respect to the classical case. On the other hand, in the J\ — J2 
Heisenberg model (Sec. 4), quantum fluctuations turn out to be strong enough to 
melt the antiferromagnetic Neel order, driving the ground state into a non-magnetic 
phase of purely quantum-mechanical nature. 

2. Spontaneous Symmetry Breaking in a Quantum Antiferromagnet 

In this section we will review the basic concepts concerning the mechanism of 
the spontaneous symmetry breaking in a quantum antiferromagnet that will be used 
in the present work to investigate the ground-state properties of the triangular and 
of the J 1 — J2 Heisenberg models. Starting from the Lieb-Mattis theorem for the 
bipartite Heisenberg antiferromagnet, we will introduce some general features of the 
finite-size spectrum of the Heisenberg Hamiltonian. In particular, we will focus on 
the importance of the structure of the low-lying excited states, explaining also how 
the finite-size ground-state properties can be consistent with a broken symmetry in 
the thermodynamic limit. 

2.1. The Lieb-Mattis property 

One of the few rigorous results on the ground-state properties of the Heisen- 
berg model on a bipartite lattice is the Lieb-Mattis theorem. Here we will repro- 
duce the demonstration of this important result, following the paper by E. Lieb 
and D. MattisEj who extended and generalized the original results obtained by W. 
MarshallB 

Let us consider the Heisenberg Hamiltonian 

(id) 

where the sum runs over all the bonds on a e?-dimensional bipartite lattice, Si are 
spin-s operators; Jy is the (symmetric) exchange matrix such that Jy > 0, if i 
and j belong to different sublattices, and Jy < 0, otherwise. We will assume that 
the Hamiltonian cannot be split into sets of noninteracting spins, restricting also 
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for simplicity, to the case in which the number of sites of the two sublattices is the 
same. In the following N will denote the total number of sites of the lattice. 

Since the Hamiltonian (||) commutes with all the three components of the total 
spin operator, 

S = ^Si, (4) 

i 

it is known from the theory of the angular momentum that we can construct two 
operators which commute with each other and with 7i. For example, choosing the 
quantization axis along the z-direction, we can consider the total spin squared, S 2 , 
and its component along the z-axis, 5 Z , whose eigenvalues 5(5+1) and M are good 
quantum numbers for the eigenstates of the Hamiltonian, 

|^ n > = KS,M,---> , (5) 

such that H\ip n ) = E n \\jj n ). If the couplings Jij are translationally or rotationally 
invariant, also the lattice momentum and the eigenvalues of the generators of the 
crystal point group label the eigenstates of the Hamiltonian. However this restric- 
tion is not needed to derive the following results. 



2.1.1. Marshall-Peierls sign rule 

In the hypothesis stated above, the first strong result one can prove is about the 
signs of the coefficients of the expansion of the ground state of (H) in the so-called 
Ising basis whose states are specified by assigning the value of the 5f at each lattice 
site, i.e., \x) — Yii=i l m «) with S 2 |m,) = s(s + l)|mj) and 5?|mj) = m,-|m,-) where 
| mi | < s. Within this basis it is easy to distinguish between the subspaces with 
different values of the projection of the total spin on the z-axis, M. In fact, in order 
to restrict to a particular M sector, one has to use only the states of the basis, |a;), 
such that J2i m i = M. In addition, sorting the basis in order to group together the 
states with the same M, the Hamiltonian matrix assumes a simpler block-diagonal 
form, i.e., a block for each M sector. Let us restrict therefore to a particular M 
subspace. 

The first step of the proof is to perform a unitary transformation, 



exp 



+ (6) 



ITT 



whose physical meaning is to flip the quantization axis on the B sublattice. This 
defines a spatially varying reference frame pointing along the local Neel direction 
[Fig. H (a)]. The transformed Hamiltonian then results 



U^lrLU = H d +H , (7) 



where the diagonal part, 



TLd — 2_. JijSfSj i (8) 

(id) 
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is invariant in the new representation, while the off-diagonal spin-flip term, 

Ho = -\Y.MStS- + h.c) , (9) 

acquires instead an overall minus sign. Therefore, in this representation, the Hamil- 
tonian has non-positive off-diagonal matrix elements. In this case, using the hy- 
pothesis that the Hamiltonian cannot be split_mto sets of noninteracting spins one 
can demonstrate (Perron-Frobenius theorems) that the ground-state expansion 
over the chosen basis \iPm) — fx\x) has non-vanishing positive amplitudes, i.e., 

U>o. 

The latter result has two important consequences: in each M subspace the 
ground state of the Hamiltonian (^) i) is non-degenerate, and ii) obeys, in the 
original representation, the well-known Marshall- Peierls sign rule, i.e., 

hM =u\M = Y, e ™ N(x) f*\ x ) = Y.(- 1 ) N(x) fx\ x ) > ( 10 ) 

X X 

where N(x) — ^2j,p^(s+mi), an( l the sum is restricted to the configurations |x) such 
that J^. m, — M. Notice that the unessential constant term s in the definition of U 
have allowed us to write a real ground-state wavefunction. In particular, for s = 1/2, 
N(x) is simply the number of spins up on the B sublattice, say Nf(x). In this case, 
the ground-state projections on two Ising configurations differing for a single spin 
flip have opposite signs. Notice that the rotational invariance of the Hamiltonian 
has never been used in the proof so that this result is valid in general even in 
presence of an easy-plane anisotropy, i.e., for the more general XXZ Hamiltonian, 
or in presence of an arbitrary magnetic field in the z direction. 

In contrast, for frustrated spin systems, like the J\ — Ji model and the Heisenberg 
triangular antiferromagnet, the same proof does not hold. In fact, in both systems 
the off-diagonal part of the Hamiltonian cannot be made non-positive defined by 
any known unitary transformation. In the J\ — J2 model, the unitary transformation 
(|J) does not change the signs of the next-nearest-neighbors spin-flip term. In the 
triangular antiferromagnet, under the transformation corresponding to the one in 
Eq. © 



exp 



(11) 



where B and C label two of the three sublattices as shown in Fig. |l]?the Hamiltonian 
reads 

U^HU = -jY.&Sj + h.c.) + jJ2S-S* 

(id) (i,j) 
+ i^Cy(^-k) (12) 



a As the unitary transformation in Eq. (|3|) this mapping defines a spatially varying coordinate 
system pointing along the local Neel direction. 
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Table 1. Weight of the states satisfying the Marshall-Peierls sign-rule, (s), in the ground state of 
the J1 — J2 Heisenberg model. Data are reported for s = 1/2, N = 32 and N = 36, and various 
values of the frustration. 



J2/J1 


0.1 


0.2 


0.3 


0.4 


0.5 


N = 32 


1 


1 


~ (1 - 10" 8 ) 


~ 0.9998 


0.973 


N = 36 


1 


1 


~ (1 - 10" 8 ) 


~ 0.9995 


0.961 



with dj = ±\/3/4. The transformed Hamiltonian then displays an extra current- 
like term which is off-diagonal and has no definite sign. Therefore for frustrated 
spin systems the Marshall-Peierls sign rule cannot be demonstrated and in general it 
does not hold exactly. In addition, the impossibility to find a unitary transformation 
allowing us to map the Hamiltonian into an operator with non-positive defined off- 
diagonal matrix elements has also dramatic consequences on the computability of 
such frustrated systems with standard quantum Monte Carlo methods. This is the 
origin of the so-called sign-problem instabilityEj'tB 

However, as originally pointed out by Richter and co-workersJEI the Marshall- 
Peierls sign-rule survives in the J\ — J2 model in very good approximation up to 
relatively large values of the frustration. By means of the Lanczos exact diagonal- 
ization technique, for s — 1/2, N — 32 and 36, we have calculated the weight of the 
states satisfying the Marshall-Peierls sign-rule in the expansion of the (normalized) 
exact ground state |^o); namely: 

(s) = \Mx)\ 2 (-l) Nr(x) sg4Mx)} , (13) 

x 

with ip(x) = (x\ip) and the notations introduced in this section 1 _Dur results, shown 
in Tab. [j], put further evidence to the previous findings of Ref. cJ and indicate that 
the Marshall-Peierls sign rule is verified almost exactly up to J2/J1 — 0.3 -r 0.4, 
even for N = 36. Moreover, even if the average sign (s) eventually vanishes in the 
thermodynamic limit, its small size dependence suggests that this property is likely 
to be conserved also for the lattice sizes (N ~ 100) presently accessible with the 
stochastic numerical techniques used in this work. A reasonable guess on the phases 
of the exact ground state is in general very useful to improve the efficiency of the 
approximated quantum Monte Carlo techniques that have to be used in presence of 
the sign problem. 0e1 

In the triangular case, instead, the phases obtained by applying the operator 
( |Tl| ) to a state with positive-definite amplitudes on the Ising basis, as in Eq. (|l0|), 
are very far from being exact, especially in the spin-isotropic limit (see Sec. 3.2.1). 

2.1.2. Ordering of the energy levels 

It follows from the spin rotational invariance of the Hamiltonian (^) that each 
energy level, E(S), belonging to the total spin S subspace is (2S , + l)-fold degenerate: 
in any S sector there is a degenerate level for each value of M in the range — S < 
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M < S. Therefore, in a given M subspace every energy eigenstate with total spin 
S > \M\ must be contained. In the hypothesis stated above, it is possible to prove 
that the lowest energy in each M subspace belongs to S = M, i.e., it has the 
minimum total spin allowed. 

In order to prove this result we will show that the ground state of H in an 
M subspace is not orthogonal to the ground state of a rotational invariant soluble 
Hamiltonian, which is known to belong to the S = M sector. Therefore, so does the 
former since two eigenfunctions having different quantum numbers are in general 
orthogonal. Let us consider the infinite-range Heisenberg Hamiltonian on a bipartite 
lattice 

n oc = j Si • §i , (14) 

ieA,j£B 

with J positive constant. This Hamiltonian is rotationally invariant and exactly 
soluble since it is equivalent to a two spin problem: 

Hoc = J(S 2 - - S%) (15) 

where S A and are the total spin squared on the A and B sublattices, respectively. 
The eigenvalues of this special Hamiltonian are 

E 00 (S) = ^[S(S+1)-S A (S A + 1)-S B (S B + 1)] , (16) 

and are monotonically increasing with the total spin S. Then, the ground state of 
Hoc, in each M subspace, has S = M total spin. Moreover, both Tt and Ttoo satisfy 
the requirements for the Marshall-Peierls sign rule. Hence, their ground state in 
any M sector are not orthogonal since their overlap involves the sum of positive 
numbers. It follows that they have necessarily the same total spin quantum number 
S = M. Therefore in a given M subspace the lowest energy of the Heisenberg 
Hamiltonian has the minimum total spin allowed. 

This implies in turn that E(S) < E(S+l). In fact, among the degenerate eigen- 
functions with E(S+l), there is a representative in the M — S subspace. The latter 
has not the minimum total spin allowed for that subspace and therefore it has an 
energy higher than E(S). This proves that the energy levels of the Heisenberg an- 
tiferromagnet (^|) increase monotonically with the total spin and, in particular, that 
the absolute ground state is a singlet and non- degenerate (Lieb-Mattis property). b 

The above proof on the ordering of the energy levels is a direct consequence of the 
Marshall-Peierls sign rule and therefore it breaks in presence of frustration. However 
the Lieb-Mattis property turns out to be verified even for frustrated spin systems. 
In particular, for symmetry reasons, the ground state on any finite size of the 
spin-isotropic Heisenberg antiferromagnet is believed to possess all the symmetries 
of the Hamiltonian and in particular to be a singlet, rotationally invariant, and 

^In fact, having S = 0, it has only a representative in the M = subspace 
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non-degenerate. cB Even if there is no rigorous theorem proving this property in 
general, the latter turns out to be true on a finite size whenever the cluster is 
large enough, it has an even number of sites, an d t he boundary conditions do 
no frustrate the antiferromagnetic long-range order.Ej In any case, however, these 
symmetry properties concern in general the ground states on finite sizes only. In the 
thermodynamic limit, the situation can change drastically if there is no gap in the 
excitation spectrum. In this case, in fact, a family of excited states collapses onto 
the ground state and may break its symmetric character. This will be illustrated 
in the following sections. 

2.2. Order parameters and susceptibilities 

A zero-temperature spontaneously broken symmetry occurs when the ground 
state has a lower degree of symmetry than the corresponding Hamiltonian. In 
this case, one can define an extensive operator, O, breaking some symmetry of 
the Hamiltonian and such that the so-called order parameter, i.e., the ground- 
state expectation value m = (ipo\0\ipo)/N , has a finite value. In general, whenever 
the symmetry-breaking operator O does not commute with the Hamiltonian, the 
symmetry breaking can happen only in the thermodynamic limit. In fact, in that 
case, the ground-state expectation value of O is zero on any finite size by symmetry. 
This will be the case for the symmetry-breaking operators considered in this paper. 

The occurrence of a spontaneously broken symmetry can be detected by adding 
to the Hamiltonian Ti. an ordering field 5: 

H s =H-6d. (17) 

Since on a finite size the ground-state expectation value of O vanishes for 6 = 0, 
the ground-state energy per site has corrections proportional to S 2 , 

e(S) ~ e - l - X o& 2 , (18) 

Xo being the (positive-definite) generalized susceptibility associated to the operator 
O, namely: 

Xo = ^(HO(E - Hy 1 d\i> ) , (19) 

where Eq is the ground-state energy of Ti.. 

If symmetry breaking occurs in the thermodynamic limit then 

Jim lim i<V>o|d|Vo) = m ^ 0, (20) 

5— *0 N—>oc JV 

and the finite-size susceptibility has to diverge with the system size. In fact, by the 
Hellmann-Feynman theorem, the ground-state expectation value of O at finite field 
is (0}s/N = —de(S)/dd, so that, if symmetry breaking occurs in the thermodynamic 
limit, an infinitesimal field 5 must give a finite (0}$/N ~ xo$ implying that the 
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susceptibility has to diverge in the thermodynamic limit. Moreover, it is possible to 
show that the finite-size susceptibility must diverge at least as the volume squared 
TV 2 . This will be proven in the next section. 

2.2.1. Exact bounds on the susceptibilities 

Since on a finite size (ifio\0\ipo) = 0, it is convenient to introduce as the order 



parameter the quantity p — y (ipo\0 2 \ipo}/N 2 . The latter is finite in general on any 
finite size and extrapolate to a finite value in the thermodynamic limit in presence 



of long-range order, i.e., whenever m, given by Eq. (20), is finite. In this section we 
will show that, whenever symmetry breaking occurs in the thermodynamic limit, 
the corresponding susceptibility must diverge as N — > oo and, in particular, it is 
bounded from below by the order parameter times the system volume squared, 
namely \o > const p A N 2 . 

Let us define the following decomposition: 

P 2 = ^2^o|0 2 |Vo) = ^El^l (5 l^)| 2 = ^ fduS(u) (21) 



with 



S(co) = ^J2\^om n )\ 2 S(cj~u Jn ) , (22) 



N 

n=£0 



where we have introduced a complete set of eigenstates of the Hamiltonian \ifj n ) with 
eigenvalues E n , we have used the symmetry of the ground state (i.e., (ipo\0\ifjo) = 0) 
and set u> n = E n — Eq. By the Cauchy-Schwartz inequality we have: 



< 



-, 1/2 

dui wS(u>) duiu>~ S(u>) . (23) 



Now, 



fdw^S{w) = i£ -L|(V> |6|^)| 2 = ^ , (24) 

where xo by Eq. ( |l~9| ) is the susceptibility associated to the operator O. In addition 
it is straightforward to show that 

[duwsw - l^c„|(v, |o|^>l 2 = ^{[d,[n,d]}) = |, (25) 

so that, using Eqs. ((||) and (pij), we have 
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Hence, we have obtained the following lower bound for the susceptibility: 

Xo > -^N 2 . (27) 
Jo 

Therefore, if the order parameter is finite in the thermodynamic limit, the suscepti- 
bility must diverge at least as the volume squared, provided fo is a constant. This 
happens whenever the commutator of TL and O is an extensive quantity as it is 
the case for the magnetization in the Heisenberg model and for all the symmetry- 
breaking operators treated in this work. 

Moreover it is also possible to construct an upper bound for the generalized 
susceptibility xo associated to a symmetry-breaking operator O. In fact, using 
Eqs. (|I|) and (|2l]), we have 

^ = ^E^I^o|0|^)| 2 <^ [duS(u,) = ??£ , (28) 

where A is the energy gap between the ground state and the first excitation. An 
energy gap in the thermodynamic excitation spectrum is therefore incompatible 
with a spontaneously broken symmetry. The physical meaning of this quite general 
result is the following: in presence of a gap in the excitation spectrum, the ground 
state, which has generally all the symmetries of the Hamiltonian on any finite size, 
has clearly no mean to develop a spontaneously broken symmetry in the thermo- 
dynamic limit. In this case, the susceptibility is bounded [Eq. (pg})]. In contrast, 
in presence of a gapless excitation spectrum, a family (or a tower) of excited states 
can collapse in the thermodynamic limit onto the ground state and can break its 
symmetric character. In fact, these states acquire in general a phase factor under 
some operation of the symmetry group of the Hamiltonian and they can give rise to 
a symmetry-broken superposition. Whenever this happens, the related susceptibil- 



ity must diverge as the volume squared [Eq. (p7h]. In this case from Eqs. (27) and 



8|) we get a remarkable relation for the size-dependence of the spin gap, namely 



A < Aat-i . (29) 
2p z 



The mechanism underlying the spontaneous symmetry breaking leading to the 
onset of long-range Neel order in the thermodynamic ground state of a quantum 
antifcrromagnet will be discussed in more detail in the following section. 

2.3. Neel order and Anderson's towers of states 

As we have already noticed, the occurrence of a spontaneous symmetry breaking 
in the thermodynamic ground state, can be evidenced from the structure of the 
finite-size energy spectrum. 

On any finite size, the ground state of a quantum antifcrromagnet is generally be- 
lieved to be a singlet, rotationally invariant and non-degenerate, i.e., non-magnetic. 
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This is rigorously stated by the Lieb-Mattis theorem only for the Heisenberg square 
antifcrromagnet but nonetheless it can be numerically verified on small clusters also 
in presence of frustration for the triangular Heisenberg antiferromagnettLS and for 
the J\ — Ji model itself.cj Therefore, as originally pointed out by Andersonjiil the 
spontaneously symmetry breaking mechanism necessarily involves the low-energy 
portion of the excitation spectrum: in particular, a whole tower of states has to 
collapse in the thermodynamic limit onto the ground state faster than the low-lying 
excitations involving a spatial modulation of the classical Neel state (the so-called 
magnons). In fact, since in general these states acquire phase factors under rota- 
tions in the spin space, they can sum up to a nontrivial state in which the spins 
point in a definite direction, giving rise to Neel-like long-range order. 

In particular, it is well knownEZhEJ that in this case the low- lying excited states 
of energy E(S) and spin S are predicted to behave as the spectrum of a free quantum 
rotator (or quantum top) as long as S <C y/N, 

E(S) -E = ^±^1 , (30) 

where Eq = E(0) is the energy of the ground state, \ipo), and I is known as the 
momentum of inertia per site and is an intensive quantity. 

This equation, which is in agreement with the bound (|29|), can be justified 
in a semi classic al picture of the long-range ordered ground state of a quantum 
antifcrromagnet .E3 To this purpose, let us consider the nearest-neighbor Heisen- 
berg antifcrromagnet on the square lattice, and separate in the Fourier transformed 
Hamiltonian the k = 0, and Q = {it, it) contributions (i.e., the only ones allowed 
by the sublattice translation invariance of the classical Neel state) from the others: 

H = H +V, (31) 

where 

7io = ^(S 2 -Si-S|) , (32) 

S 2 is the total spin square, and S\ and are the total spin square of the A and 
B sublattices, respectively; 

V = 2J J2 7kS k -S_ k , (33) 

k#0,Q 

where Sk = l/V^/V^i ^ exp(k • rj), is the position of the site i, and 7k = 
(cos k x + cos k y )/2. For each value of the total spin S the lowest eigenstate of Hq is 
the classical-like state fully polarized on each magnetic sublattice with energy: 

* W __» (W + 4) + «»£±i>. ,34, 



This is in agreement with the quantum top law (|30|). Of course, since and Sfj 
do not commute with V, such eigenstates are not eigenstates of 7i. Nonetheless we 
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can look at them as the first approximation to the low-lying excited states in each 
5 sector. The perturbation V dresses these classical-like states with quantum fluc- 
tuations decreasing the average value of the sublattice magnetization and lowering 
their energy towards the exact result. However, as long as the energy scale of these 
states is well separated by the low-lying excitations with k 0, Q such renormal- 
izations are not expected to modify the behavior (|3^) . These excitations are known 
as magnons or spin-waves (see Sec. 3.1), and involve a spatial modulation of the 
classical Neel state. In the Heisenberg antiferromagnet the dispersion relation of 
the softest magnons is linear in the wavevector so that, in two dimensions, they 
have an energy scaling as 1/ yN. This implies that the constraint on the value of 
5 for the validity of Eq. © is 5 < VN. 

The quantum top law is similar to the definition of the uniform spin susceptibility.EU 
The latter, in fact, can be calculated, by taking first the infinite- volume limit of the 
energy per site e(m) = E(S)/N at fixed magnetization m = S/N and then letting 
m — > in the expansion 

o 

in 

e(m) = e + — , (35) 
2x 

which is quite similar to Eq. (|30|). However an identification between / and % 
is possible if the excitation spectrum smoothly connects the low-energy portion 
which corresponds to total spin 5 pi 0(1), with the regime of macroscopic spin 
excitations: 5 ~ mN (with m <C 1).E3 This is an highly nontrivial statement which 
is actually verified by the underlying low-energy effect hzB_rm dfiLof the quantum 
antiferromagnet, known as nonlinear a model (NLcrM).^!^'^^!^ Therefore, the 
quantity 

Xs 5(5 + 1) 

must approach the inverse of the spin susceptibility for infinite size and for any spin 
excitation with 5 <C N. 

2.4. Resonating Valence Bond states 

A simple and clear picture of a non-magnetic ground state canine given in terms 
of the so-called Resonating Valence Bond (RVB) wavefunctions £.3 Here, for sim- 
plicity, we will restrict ourselves to the case of the spin-half square antiferromagnet 
even if these states can be used also for a generic value of the spin s (Ref. £3) and 
different lattice geometries. E3 

The RVB wavefunctions are linear superpositions of valence bond states in which 
each spin forms a singlet bond with another spin on the opposite sublattice. These 
states form in general a (overcomplete) basis of the 5 = subspace so that any 
singlet wavefunction can be represented in terms of them. In particular, with such 
RVB wavefunctions it is possible to describe either a long-range ordered or a non- 
magnetic state by varying the bond-length distribution.LJ In order to clarify this 
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Figure 3: An example of a spin liquid (a) and of a symmetry-broken (b) non- 
magnetic RVB state. Each stick represents a singlet bond. 

point, let us consider the following class of RVB wavefunctions for a system of N 
spins: 

\iPrvb) = ^2 H r i)---h(r n ) (i 1 j 1 )...(i n ,j n ) , (37) 

where n — N/2, r rn is the distance between the spins forming the m th singlet 
bond (i m jm), and h(r m ) is a bond weight factor, function of its length. The latter 
wavefunction has no long-range order whenever the short-ranged bonds are the 
dominant one in the superposition (^). More precisely, it has been numerically 
shown by Liang, Doucot and Andersonlll that the RVB state (^) has no long- 
range antiferromagnetic order for bonds that decay as rapidly as h(r) ~ r~ p , with 
p > 5. Instead, if the weight functions decay slowly enough with the length of the 
bond, then the RVB wavefunction has a finite value of the thermodynamic order 
parameter squared. In particular, if the weight factors h(r) are independent of the 
bond length, the RVB wavefunction is the projection of the Neel state onto the 
singlet subspace. 

Therefore, the simplest physical picture of a non-magnetic ground state can 
be given in terms of a RVB wavefunction with short-ranged bonds. In addition, 
such bonds can be either homogeneously spatially distributed on the lattice, with 
short-range correlations among each other [spin liquid) [Fig. |^ (a)], or they can 
break some symmetries of the Hamiltonian, with the dimers frozen in some special 
patterns [Fig. | (b)]. In Sec. 4.3 we will provide a possible example of the latter 
situation. 

3. The Triangular Heisenberg Antiferromagnet 

Historically the antiferromagnetic spin-1/2 Heisenberg model on the triangular 
lattice was the first pr op osed Hamiltonian for a microscopic realization of a non- 
magnetic ground state.oE3 This is due to the fact that in this system the usual 
antiferromagnetic alignment between spins is hindered by the geometry of the lattice 
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so that the minimum energy configuration, the 27r/3 Neel state (Fig. |l|), has an 
energy twice larger than that on the square lattice and therefore is far less stable. 

However, despite a strong theoretical effort, the question of whether the com- 
bined effect of frustration and quantum fluctuations in the ground state of the 
triangular Heisenberg antiferromagnet favors a disordered RVB state or long-range 
Neel type order has been under debate for many years.0ElBil~il 

Spin-wave calculations0il predict an important reduction (by about one-half) 
of the siiblattice magnetization by quantum fluctuations. In addition, perturbation 
theory,E3 series expansions ,c£l and high-temperature calculationscll suggest that the 
spin-wave calculations possibly underestimate the renormalization of the order pa- 
rameter, but do not come to a definite conclusion about the nature of grc 
From the numerical point of view, exact diagonalization (ED) results,!! 
which are limited to small lattice sizes, have been interpreted both againstEil and 
in favorlij of the presence of long-range Neel order in the thermodynamic ground 
state. In any case, this approach, combined with the careful analysis of the symme- 
try properties of the low-energy excited states proposed by Bernu and co-workers ,E3 
have provided very important evidence pointing towards a magnetic ground state: 
the spectra of the lowest energy levels order with increasing total spin, a remi- 
niscence of the Lieb-Mattis theorem (see Sec. 2.1) for bipartite lattices, and are 
consistent with the symmetry of the classical order parameter. However, these very 
clear ED results cannot rule out that for large sizes quantum fluctuations could 
drive the system into a non-magnetic phase and therefore cannot be considered as 
conclusive. In addition, standard stochastic numerical methods, which usually allow 
one to handle large samples, clash with the sign problem numerical instability so 
that a definite answer on the ground-state properties of the triangular Heisenberg 
antiferromagnet has been lacking for many years. 

In this section we will tackle the problem of the existence of long-range Neel order 
in the ground state qfthe triangular Heisenberg antiferromagnet using the finite- 
size spin-wave theory,E2l ED, and Green function Monte Carlo (GFMC) methods. In 
the first part we apply the finite-size spin- wave theory to the triangular Heisenberg 
antiferromagnet, we then show how to construct within this framework the low- 
lying excited states, and finally derive a simple spin-wave variational wavefunction. 
The good agreement between the ED results and the finite-size spin- wave theory will 
support the reliability of the spin- wave expansion in describing not only the ground- 
state properties but also the low-energy spin excitations of the Heisenberg model 
even in presence of frustration. The second part will deal with the quantum Monte 
Carlo results. Data for the spin gap and for the antifcrromagnetic order parameter 
will be presented for fairly large system sizes (up to 144 sites), providing a robust 
evidence for a gapless excitation spectrum and for the existence of long-range Neel 
order in the thermodynamic ground state of the model. 

3.1. Finite-size spin-wave theory 

Several attempts to generalize spin- wave theory to finite sizes can be found in the 
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literaiure.EJi£j'E2l Here we will follow the method proposed by Zhong and Sorella in 
Ref. E3 which allows one to deal with finite clusters avoiding the spurious Goldstone 
modes divergences in a straightforward way, and, in particular, witho ut imposing 
any ad hoc holonomic constraint on the sublattice magnetization.1 



15SU56 



3.1.1. Application to the triangular antiferromagnet 

Assuming the classical Q = (47r/3,0) magnetic structure lying in the xy plane, 
the first step of the derivation is to apply the unitary transformation given by 



Eq. (11), which defines a spatially varying coordinate system (x'y'z') in such a way 
that the x'-axis points on each site along the local Neel direction. The transformed 
Heisenberg Hamiltonian reads: 

U^HU = J [cos (Q • r id ) (sf sf + Sf'sf) 



S m(Q-r id )(Sfsy'-Sy'Sf) + SfSf (3M 



where J is the (positive) exchange constant between nearest neighbors, the indices 
i,j label the points rv and rj on the iV-site triangular lattice, r^.j = — r 3 , and 
the quantum spin operators satisfy |S;| 2 = s(s + 1). In the new reference frame 
the spins in the classical configuration are ferromagnetically aligned so that, using 
Holstein-Primakoff transformation for spin operators to order 1/s, 



being a and a' the canonical creation and destruction Bose operators, after some 
algebra the Fourier transformed Hamiltonian results: 



3Js ^ [A k a k a k + ^Bk(a k a-k + akfl-k) (40) 



k 

where E cl = —3Js 2 N/2 is the classical ground-state energy, 

^k = l+7k/2 , 5 k = -3 7 k/2 , (41) 

7k=2 [cos(k x ) + 2 cos {k x /2) cos (V3 k y /2)] /z, k is a vector varying in the first 
Brillouin zone of the lattice, and z = 6 is the coordination number. The Hamiltonian 
Hsw, can be diagonalized for k ^ 0,±Q introducing the well-known Bogoliubov 
transformation, a k = ii k a! k + u k al_ k , with 

a - \ 1/2 / , \ 1/2 

^k + £k \ , „ > / A k - e k x 



where e k = U A 2 . — is the spin-wave dispersion relation. This diagonalization 
leads to: 

3Js 3Js 
Ksw = E ol + ^-J2 ( £k _ Ak ) + — 12 e k(a k «k + al k a-k )• (43) 

k#0,±Q k^0,±Q 
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The Goldstone modes at k = and k = ±Q instead are singular, and cannot 
be diagonalized with a Bogoliubov transformation. For infinite systems such modes 
do not contribute to the integrals in Eq. (^), but in the finite-size case they are 
important and they must be treated separately. By defining the following Hermitian 
operators 



Qx 



i(% -a ) , 



(44) 



such that, [QaiQp] — and [Q Q ,7isw] = for a,(3 = x,y,z, the contribution of 
the singular modes, Hsm, in Eq. ( |40| ) can be expressed in the form 

A 



Hsm = —3JsA + 3Js- 



Qi + Qt 



(45) 



Then, taking into account the fact that to the leading order in 1/s, Q a — S a y / z/Ns, 
where S a are the components of the total spin, TLsm may be also rewritten in the 
more physical form 



n 



SM 



-3JsA 



A 
3J W 



(S x ) 2 + (S y ) 2 + (S z ) 2 



(46) 



which clearly favors a singlet (S 2 = 0) ground state (for an even number of sites) 
being Ao positive definite. This result is a reminiscence of the Lieb-Mattis property 
(see Sec. 2.1) which has not been demonstrated for non-bipartite lattices. Actually, a 
similar result is obtained by solving exactly the three Fourier components k = 0, ±Q 
of the Heisenberg model;E3 however, our treatment allows us to construct a formal 
expression for the spin-wave ground state on finite triangular lattices which keeps 
the correct singlet behavior. In fact, starting from the usual spin- wave ground state, 
composed by the 2ir/3 classical Neel order plus the zero point quantum fluctuations 
(i.e, zero Bogoliubov quasiparticles), 



io>= n «i 

k#0,±Q 



exp 



2 Uk 



\F) 



(47) 



with \F) = Y\i \S£ = s )> the corresponding singlet wavefunction is obtained by 
projecting |0) onto the subspace S = 0: 



oo />oo 



|V>SW>=/ da d0 e iaQ ' +i ^ Q y +i ^\0) 



(48) 



' — oo «/ — oo <J — oo 
at _i_ l at at * 



and reads IV'Sw) ~ e< * a Q a -Q +2a ° a °- l |0) . In particular the singular modes have no 
contribution to the ground-state energy which reads 



^sw — E c i H — — ^(fk — 1) , 



(49) 



20 Quantum Effects and Broken Symmetries in Frustrated Antiferromagnets 



while the computation of the order parameter requires their remotion: 



,t _ 



<(<?f') 2 )= S - 



k#0,±Q 



(50) 



For s — 1/2, the previous spin- wave calculation predicts a very good quantitative 
agreement with exact results on small clusters (N < 36) of both ground-state energy 
and sublattice magnetization.E3 The agreement is even more remarkable as far as the 
low-lying excited states are concerned, as it will be shown in the following section. 



3.1.2. Low-energy spin-wave spectrum 

In this section, we show how to construct the low-lying energy spectra E(S) 
for finite systems where S represents the total spin. So far, we have performed a 
standard spin-wave expansion whose relevant quantum number is s. Thus, com- 
puting E(S) is not straightforward and require a little more involved calculation. 
Following Lavalle, Sorella and Parola@ a magnetic field in the z-direction is added 
to stabilize the desired total spin excitation S: 



Ti-sw 



(51) 



Classically, for magnetic fields not large enough to induce a spin-flop transition, the 
new solution is the 27r/3 Neel order canted by an angle 6 along the direction of the 
field h. In order to develop a spin-wave calculation, a new rotation around y'-axis 
is performed on the spin operators and it can be proven that 7ig W takes the same 
form of Eq. (Boh with renormalized coefficients ^4k and B^: 



A* = l 



7k 



1 3 2h 2 

2 ~ Vl&j' 



, 2h 

Z6J 



(52) 



being 2h/z3J — sin0. For h ^ the only singular mode is k = (associated to the 
rotation invariance in the xy plane) and its contribution is given by 



Hsm = ~A^ + 3J^ (S z - Ns sin^) 2 , 



(53) 



which now favors a value of S z consistent with the applied field, at the classical 
level. 

The Hellmann-Feynman theorem relates the total spin S = N(S*} of the exci- 
tation to the magnetic field h as it follows: 



w—kmw-ib 




where 



E(h) 



Ecl - i( 8h) tej " 3js ~ + — E £ k- 



(54) 



(55) 
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Figure 4: Spin-wave (full dots and continuous line) and exact (empty dots and 
dashed lines) low-energy spectra as a function of |S 2 | = 5(5 + 1) for N = 
12,27,36,144 and s = 1/2. 



is the spin-wave energy in presence of the field and e£ = y (A£) 2 — (B^) 2 . In 
particular, the first term in (sh) 2 in Eq. ( |55| ) gives the classical uniform spin sus- 
ceptibility Xci = 1/9 J, while taking the whole expression the known spin- wave 
resuhrJ Xswjxd = 1 — 0.449/2s is recovered. Finally, as suggested by Lavalle and 
co-workerSjt 2 ] given the value 5, the corresponding values of h and of E{h) can be 
found with Eqs. (pi] ) and ((5^), and the energy of the spin excitation E(S) can be 
calculated, at fixed s, with a Legendre transformation E{S) — E(h) + hsSz3 

As explained in Sec. 2.3, the occurrence of a symmetry breaking in the ground 
state for N — > oo can be evidenced from the structure of the finite-size energy spec- 
tra. In particular, when long-range order is present in the thermodynamic limit, the 
low-lying excited states of energy E(S) and spin 5 are predicted to behave as the 
spectrum of a free quantum rotator (|3(]) as long as S <C V^V- Actually, on the trian- 
gular lattice the quantum-top effectiveJIamiltonian displays a correction due to the 
anisotropy of the susceptibility tensor J13 However, in the following we will consider 
only the leading contribution ( |3C| ) which depends on the perpendicular susceptibil- 
ity. Fig. |^ shows E(S) vs 5(5+1) calculated within the spin- wave theory compared 
with the exact diagonalization results of Bernu and co-workers lid Remarkably the 
spin-wave theory turns out to be accurate in reproducing the low-energy spectrum 
in the whole range of sizes. Furthermore, we can extend our calculation to the 
thermodynamic limit and observe easily the collapse of a macroscopic number of 
states with different 5 to the ground state as N — > oo. This clearly gives rise to a 
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Figure 5: Size dependence of 1/2xl (a) and of 1/2xl/2~1/2xl (b) obtained accord- 
ing to Eq. ( |56| ) using the (s = 1/2) spin- wave excitation spectra. The continuous 
line is a quadratic fit for L < 18 in (a) and a guide for the eye in (b). 

broken SU(2) symmetry ground state, as expected within the spin- wave framework. 

In addition, as explained in Sec. 2.3, whenever the quantum top law ( ^p| ) is 
verified, the quantity 

pxsP 1 =NE(S)[S(S+1)}- 1 , (56) 

should approach the physical inverse susceptibility 1 /2xsw f° r infinite size and for 
any spin excitation S <C N. This feature is clearly present in the spin-wave theory 
and it is shown in Fig. || (a) where the l/2xs is plotted for S = L = \/~N and 
approaches the predicted value (l/2xsw = 8.167), even if the correct asymptotic 
scaling 1/2\l — l/2xsw + a /L + b/L 2 turns out to be satisfied only for very large 
sizes [L > 36). Such feature is also shared by the Heisenberg antiferromagnet on 
the square lattice where a similar spin-wave analysis has allowed the authors of 
Ref. E3 to account for the anomalous finite-size spectrum resulting from an accurate 
quantum Monte Carlo calculation. Furthermore, similarly to the latter case, a non- 
monotonic behavior of 1/2xl/2 — 1/2xl [Fig.|| (b)], which should extrapolate to as 
1 / L according to the quantum top law, persists also in presence of the frustration 
within the spin-wave approximation and is likely to be a genuine feature of the 
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Heisenberg model. 

3.1.3. Spin-wave variational wave junctions 

In a GFMC calculation it is in general important to start from an accurate 
variational guess of the ground-state wavefunction. So far, many wavcfunctions 
have been proposed in the literatur and the lowest ground-state energy 

estimation was obtained with the long-range ordered type.ESc3 

The simplest starting point for constructing a long-range ordered wavefunction 
is of course the classical Neel state. Since on finite-size the ground state is expected 
to be rotationally invariant, the Neel state should be projected onto the S 2 = 
subspace. However it is in general numerically very difficult to perform the projec- 
tion onto a total spin subspace so that only the projection onto the subspace with 
S z = (a quantum number of the states of the chosen basis) is usually performed. 
Quantum correction to this classical wavefunction can be introduced by means of a 
Jastrow factor containing all the two-spin correlations 

|^)=Poexp(|^ V (i-i)^)|iV) , (57) 

where Pq is the projector onto the S z = subspace. Starting from the spin-wave 
ground state ( [48| ) it is possible to deriveS a simple variational wavefunction which 
is both accurate and easily computable also when used for importance sampling in 
a quantum Monte Carlo calculation.^ Such wavefunction is defined for any s in the 
correct Hilbert space of the spin operators and reduces for s — ► oo to the spin wave 
form @. 

To this purpose let us consider the following variational wavefunction 



\$ G ) = P cxp [1 J2 flk&$ik] \F) , (58) 



- s 

where = N^ 1 / 2 £V e~ lk r i S?. The spin- wave s — > oo limit of the wavefunction 
( |58|) can be easily carried out by means of a Hubbard- Stratonovich transformationES 
and leads to (neglecting an unessential normalization) 

|V3 G > - Po exp [ - W T^«k«- k ] (59) 

By requiring that |i/'g) reduces to the spin- wave wavefunction ( p8| ) for s — ► oo one 
obtains for gk 



«k 1 /l + 27 k 

.9k = = 1 - \ (60) 

Vk - Uk V 1 - 7k 

which is singular only for k = 0. This analysis, for the more general XXZ Hamilto- 
nian with an exchange easy-plane anisotropy apQ gives 



9k = l-\j : _ 7k ■ (61) 
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In the original (unrotated) reference frame, the Neel state \N) can be written in 
terms of \F) — Y\ i \Sf = s) by applying the inverse of the unitary transformation 

\n) =u\f) = y,u{x)\ X ) = E cx p[+^r(E^-E 5 * ) 1 1*) . ( 62 ) 



where \x) is an Ising spin configuration specified by assigning the value of Sf (or 
equivalently of Sf ) for each site, and U(x) = (x\U\x). Then, introducing a varia- 
tional parameter r\ scaling the latter potential, for s = 1/2 and in the original spin 
representation, \ipG) assumes the very simple form of Eq. (|57|), i.e., 

IV^g) = U$g) -E W W ex P [|E u ^'-i)^i]k) . (63) 

where u(r) = 1/Wj^ ^ e_ * q *5q an d the summation is restricted only to the Ising 
configurations with ^ 5,f = to enforce the projection onto the S z = subspace. 
In contrast to the linear spin- wave ground state (ff8|), which does not satisfy the 
constraint {a\dn) < 2s, the present variational wavefunction is defined in the correct 
Hilbert space of the spin operators. 

3.2. Quantum Monte Carlo calculation 

3.2.1. From Marshal- Peierls to Huse-Elser sign rule 

According to the Marshall theorem (see Sec. 2.1), for the Heisenberg antiferro- 
magnet on the square lattice, as well as on any other bipartite lattice, the classical 
part of the wavefunction by itself determines exactly the phases of the ground state 
in the chosen basis. For the triangular case, instead, the exact phases are unknown 
and the classical part is not enough to fix them correctly. In particular,^tarting from 
a long-range ordered wavefunction of the form (|57|), Huse and Elseie3 introduced 
important three-spin correlation factors in the wavefunction 

f = exp(z/3 E VjkSffyH) - (64) 

defined by the coefficients jijk — 0, ±1, appropriately chosen so as to preserve the 
symmetries of the classical Neel state, and by an overall factor f3. In particular 
the sum in Eq. ( |(34]) runs over all distinct triplets of sites i, j, k where both i and k 
are nearest neighbors of j, and i and k are next-nearest neighbors to one another. 
The sign factor 7^ = 7^ = ±1 is invariant under rigid translations and rotations 
in real space by an angle of 27r/3 of the three-spin cluster i,j,k, but changes sign 
under rotations by 7r/3 or 7r. The resulting wavefunction reads therefore: 

|V> G ) = E 6x p (I E < l ~ fi s t s i) l*> > ( 65 ) 
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Figure 6: Average sign, overlap square and ground-state energy per site obtained 
for N — 36 using the variational wavefunction of Eq. (^5|), with (full dots) and 
without (empty dots) the triplet term of Eq. (Q), as a function of the easy-plane 
anisotropy a. The calculations were performed by summing exactly over all the 
configurations and the dotted line connects the exact results. 



with a phase factor given by 

£l(x) = T(x) exp 



2-ki 

IT 



ieB 



(66) 



where T{x) — (x\T\x). Finally, since the Hamiltonian is real, a better variational 
wavefunction on a finite size is obtained by taking the real part of Eq. (|6^) . 

Although the three-body correlations of Eq. ( |64| ) do not provide the exact an- 
swer, they allow us to adjust the signs of the wavefunction in a nontrivial way 
without changing the underlying classical Neel order. In this respect it is useful 
to define an average sign of the variational wavefunction relative to the normalized 
exact ground state \^ Q ) as 



(s) = \Mx)\ 2 sgn[ip G (x)^ (x)] , 



(67) 



with ip(x) = (x\ip). We have compared the variational calculation with the ex- 
act ground state obtained by ED on the N = 36 cluster. For completeness we 
have considered the more general XXZ Hamiltonian with the exchange easy-plane 
anisotropy a, ranging from the XY case (a = 0) to the standard spin-isotropic case 
(a = 1). As shown in Fig. |(], the introduction of the three-body correlations of 
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Figure 7: Variational estimate (VMC) and mixed averages (FN, GFMCSR) of the 
ground-state energy per site E /JN and of the total spin square (S 2 ) for N = 36. 
GFMCSR data are obtained using the short-range correlation functions generated 
by H (p = 2), and H 2 (p = 7) as explained in the text. 



Eq. (|64|), although not providing the exact answer, improves the overlap square of 
the variational wavefunction on the true ground state and the accuracy of the varia- 
tional estimate of the ground-state energy as well. In particular the average sign (s) 
is very much improved by the triplet term, particularly in the spin-isotropic limit 
a — ► 1. This is crucial when the variational wavefunction is used for importance 
sampling within the modifications of the GFMC technique developed to handle the 
sign problem instability, like the FN and GFMCSR techniques. In the following the 
wavefunction ( |65|) will be used as the guiding wavefunction in our quantum Monte 
Carlo calculations. 

3.2.2. The reconfiguration scheme 

In order to study the ground-state properties we have used the GFMCSR quan- 
tum Monte Carlo technique, which allows one to release the FN approximation, in 
an approximate but controlled way. This systematic improvement introduced by 
the GFMCSR on the accuracy of the ground-state properties, is illustrated in Fig. fj], 
where we display a comparison between the estimates of the ground-state energy 
per site and of total spin square, for the N — 36, obtained with the stochastic sam- 
pling of the variationalwayefunction (|65|), the FN and the GFMCSR techniques. 
As explained in Refs. EZH13, in the appropriate limit of large number of walkers 
and high frequency of SR, the residual bias introduced by the GFMCSR depends 
only on the number p of operators used to constrain the GFMC Markov process, 
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Figure 8: Short range spin correlation functions generated by H (a,b) and H 2 
(c-g). 

allowing simulations without numerical instabilities. In principle the exact answer 
can be obtained, within statistical errors, provided p equals the huge Hilbert space 
dimension. In practice it is necessary to work with small p and an accurate selection 
of physically relevant operators is therefore crucial. As can be easily expected, the 
short-range correlation functions SfSj and (S^ S^+S^ Sj) contained in the Hamil- 
tonian (p — 2) give a sizable improvement of the FN ground-state energy when they 
are included in the SR procedure. In order to be systematic, we have included in 
the SR also the short-range correlations generated by n 2 (p = 7), averaged over all 
spatial symmetries commuting with the Hamiltonian. These local correlations (see 
Fig. D) are particularly important to obtain quite accurate and reliable estimates 
not only of the ground-state energy but also of the mixed averageEZrtLS of the to- 
tal spin square S 2 . In particular it is interesting that, starting from a variational 
wavefunction with no definite spin, the spin rotational invariance of the finite-size 
ground state is systematically recovered by means of the GFMCSR technique (see 
lower panel of Fig. [?]). 

Having obtained an estimate for the ground-state energy, at least an order of 
magnitude more accurate than our best variational guess, it appears possible to 
obtain physical features, such as a gap in the spin spectrum, that are not present 
at the variational level. For instance in the frustrated J1 — J2 Heisenberg model 
(see Sec. 4.2), with the same technique and a similar accuracy, a gap in the spin 
spectrum is found in the thermodynamic limit, starting with a similar ordered and 
therefore gapless variational wavefunction. 

3.2.3. Ground-state energy and spin gap 

In presence of Neel long-range order, being the magnon dispersion relation linear 
in the wavevector k. the leading finite-size correction to the ground-state energy per 
site is (D(N~ 3 / 2 )x3 this is clearly shown by the behavior of the finite-size spin- wave 
results in Fig. ffl. In the same figure the size scaling of the estimates of the ground- 
state energy per site obtained with the VMC, the FN and the GFMCSR (p = 7) 
techniques is also reported. The predicted size scaling, fulfilled of course by the 
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Figure 9: Ground-state energy per site eo = Eq/N, in unit of J, as a function of the 
system size, obtained with VMC (full triangles), FN (empty dots) and GFMCSR 
with p — 7 (full dots) techniques. Spin-wave size scaling is assumed and short- 
dashed lines are linear fits against 1/iV 3 / 2 . The long-dashed line is the linear spin- 
wave prediction, the emptv triangle is the N — 36 ED result and the empty squares 
are data taken from Ref. 

variational wavefunction (|6^), is also preserved within the FN and the more accurate 
GFMCSR technique, thus providing a first clue on the existence of long-range Neel 
order in the thermodynamic ground state of the model. The quality of wr results 
is similar to the variational one obtained by Sindzingre and co-workersJ23 using a 
long-range ordered RVB wavefunction (Sec. 2.4). The latter approach is almost 
exact for small lattices, but the sign problem is already present at the variational 
level, and the calculation has not been extended to high statistical accuracy or to 
N > 48. Our best estimate is that in the thermodynamic limit the ground-state 
energy per site is eo = —0.5458 ± 0.0001 in unit of the exchange coupling. 

In the isotropic triangular antiferromagnet, the gap to the first spin excitation 
is rather small. Furthermore, for the particular choice of the guiding wavefunction 
(|65|), the translational symmetry of the Hamiltonian is preserved only if projected 
onto subspaces with total S z multiple of three. Then, we have studied the gap 
to the spin 5 = 3 excitation as a function of the system size. Technically, within 
our numerical framework, such a spin gap can be evaluated by performing two 
simulations in the S z — and S z — 3 subspaces. This can be easily done by 
restricting the sampling to the configurations |x) in Eq. ([35]) with the desired value 
of S z . In this case the potential v(r) used was the same in both subspaces and the 
variational parameter rj was found by optimizing the energy in the S z — subspace. 

As it is shown in Fig. [l^, for the lattice sizes for which a comparison with ED 
data is possible, the spin gap estimated with the GFMCSR technique is nearly 
exact. The importance of extending the numerical investigation to clusters large 
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Figure 10: Size scaling of the spin gap to the S 1 = 3 excitation obtained with 
the GFMCSR (p = 7) technique (full dots). Empty triangles are the ED results, 
the long-dashed line is the linear spin-wave prediction, the dotted line is the linear 
extrapolation of the N — 12, 36 exact results and the solid line is the least-squares 
fit of the GFMCSR data for TV > 36. 

enough to allow a more reliable extrapolation is particularly evident in the same 
figure, in which the TV = 12 and 36 exact data extrapolate linearly to a large 
finite value. This behavior, is certainly a finite-size effect and it is corrected by 
the GFMCSR data for N > 48, suggesting, strongly, a gapless excitation spectrum 
[(E 3 -E )/J = 0.002 ±0.01]. 

3.2.4. Staggered magnetization 

The GFMC allows us to obtain a very high statistical accuracy on the ground- 
state energy, but does not allow us to compute directly ground-state expectation 
values {tp \0\ipo) . A straightforward way to calculate such expectation values is to 
use the Hellmann-Feynman theorem. In fact, if the Hamiltonian is perturbed with 
a term —AO the first order correction to the ground-state energy is, by standard 
perturbation theory, 

E(X) = E - \(-4>o\0\^o) ■ (68) 

As a consequence it is possible to evaluate the ground-state expectation value 
(ipo\0\ipo) — —dE(X)/d\\\=o estimating the limit 

(V>oWo) = -lim E(A) ~ jE ° (69) 

A — >0 A 

with few computations at different small A's. 

A further complication for nonexact calculations like the FN or GFMCSR, is 
that if the off-diagonal matrix elements (x'\0\x) of the operator O (in the chosen 
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Figure 11: Size scaling of the order parameter: VMC (full triangles), FN (empty 
dots), GFMGSR (full dots), exact data (empty triangles) and finite-size linear spin 
wave theorjtD (empty squares). The inset displays the A — > extrapolation for 
N > 12. Lines are quadratic fits in all the plots. 

basis) have signs which are opposite to those of the product s/'gC^OV'gOe)) they 
cannot be handled exactly within the FN because the addition of such a pertur- 
bation to the Hamiltonian changes the nodal surface of the guiding wavefunction. 
In that case, in fact, the effective FN Hamiltonian associated to the unperturbed 
Hamiltonian is also affected by the presence of the field and this leads naturally 
to the breakdown of Eq. (p8|). A way out of this difficulty if to split the operator 
O into three contributions: O = D + + + 0~, where + (0~) is the operator 
with the same off-diagonal matrix elements of O when they have the same (oppo- 
site) signs of ipaix'^dx), and zero otherwise, whereas D is the diagonal part of 
O. Then, we can add to the Hamiltonian a contribution that does not change the 
nodes: H(\) = H - \{D + 2 6+) for A > and H(\) = H - X(D + 2 6~) for A < 0. 
Then the expectation value of the operator O can be written as 

(«„)=limtt^. (70) 

With this method, using the FN and GFMCSR techniques, we have calculated 
the order parameter 



(71) 
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where Ai 2 is the sublattice magnetization squared. Our results are plotted in 
Fig. [nj. For the order parameter the inclusion of many short-range correlations in 
the SR is not very important. Then, in order to minimize the numerical effort, we 
have chosen to put in the SR conditions the first four correlation functions shown 
in Fig. ^ the order parameter itself and S 2 . While the FN data extrapolate to a 
value not much lower than the variational result, the GFMCSR calculation provides 
a much more reliable estimate of the order parameter with no apparent loss of 
accuracy with increasing sizes. In this way we obtain for m) a value well below 
the linear and the second order (which has actually a positive correction^ 2 ]) spin- 
wave predictions. Our best estimate is that in the thermodynamic limit the order 
parameter rw = 0.41 ±0.02 is reduced by about 59% from its classical value. This is 
partially in agreement with the conclusions of the finite-temperature calculationscJ 
suggesting a ground state with a_j3mall but nonzero long-range antiferromagnetic 
order and with series expansions^ indicating the triangular antiferromagnet to be 
likely ordered but close to a critical point. In our simulation, however, which to our 
knowledge represents a first attempt to perform a systematic size-scaling analysis 
of the order parameter, the value of rrv remains sizable and finite, consistent with 
a gapless spectrum. 



4. The Ji J 2 Heisenberg Model 

In the triangular Heisenberg antiferromagnet, frustration is induced by the ge- 
ometry of the lattice. The other possible origin of the frustration comes from 
competing interactions as in the so-called J1 — J2 Heisenberg model on the square 
lattice defined in Eq. (||) of the Introduction. 

In th^last few yeara^saveral studies - includin g ex act diagonalizationsjaLamall 
clusters Bj'E^I spin- waveE3E2ll£3 and Schwinger-bosonE^I calculations, seriesE3'E3'E3 and 



large- iVEJ expansions - have provided some evidence for the absence of Neel order 
in the ground state of the spin-1/2 J1 — J2 Heisenberg model for 0.38 < J2/J1 < 0.6. 
However, due to the difficulties related to the sign problem, a systematic size-scaling 
of the spin gap has been lacking for many years, and no definite conclusion on the 
nature of the non-magnetic phase has been drawn yet. In particular, an open 
question is whether the ground state_m the quantum disordered phase is a RVB 
spin liquid with no broken symmetryJ23 or if it breaks some crystal symmetries by 
dimcrizing in some special pattern (Sec. 2.4). 

In this section, we will address the problem of the quantum phase-transition to a 
non-magnetic ground state driven by frustration in the spin-1/2 J\— J2 Heisenberg 
model by means of the finite-size spin-wave theory, ED and GFMC techniques. 



4.1. Finite-size spin-wave results 

The finite-size spin-wave theory for the spin-s J\—J% Heisenberg model can be 
derived along the same lines followed for the triangular Heisenberg antiferromagnet 
in Sec. 3.1. In this section we will only show the main results for J2/J1 < 0.5, 
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assuming the two sublattice classical Neel order [Fig. |^ (a)] in the xy-pl&ne. 

Applying the unitary transformation (^) which rotates the spin quantization 
axis by a angle ir about the z-axis on one of the two sublattices, setting the or- 
der parameter along the s-axis, and using the Holstein-Primakoff transformation 
for spin operators at the leading order in 1/s, the Fourier transformed spin- wave 
Hamiltonian results as in Eq. (^fj) with E d = -zJs 2 N(l - /3)/2, z = 4, (3 = J2/J1, 

A k = 1 + p(6 k - 1) , B k = - 7k , (72) 

7k=(cos k x + cos fcy)/2, and <5 k = cos k x cos k y . 

Similarly to the triangular lattice case, the singular Goldstone modes [k = 
and k = (71", it)] cannot be diagonalized by means of the Bogoliubov transformation 
but can be recombined to give the total spin squared S 2 at the leading order in 1/s: 

Hsm = -JiszAo + Jiz^ [(S x ) 2 + (Sy) 2 + (S z ) 2 ] ■ (73) 

As seen in Sec. 3.1 this term, being Aq positive definite, favors a singlet ground 
state and implies the Lieb-Mattis property, which has been demonstrated only for 
bipartite Hamiltonians, but nonetheless can be verified numerically on finite sizes 
for the J1 — J2 Heisenberg model. 

As in the triangular case, the above analysis allows one to derive a variational 
wavefunction which is both accurate and easily computable in a quantum Monte 
Carlo algorithm. Here we will not repeat the derivation, which follows very closely 
the one for the triangular antiferromagnet, and leads to the following result for 
s = 1/2: 

x i,j 

where v(r) — l/^E q ^o e_lk r J k w ^ n 



u k - Wk V 1 - P(l - S k ) - 7k 

I x) is an Ising spin configuration specified by assigning the value of Sf for each site 
and S M {x) = (-1)^0) is the Marshall sign (Sec. 2.1), depending on the number 
N^(x) of spin up on one of the two sublattices. The summation is restricted only to 
the Ising configurations with Sf = in order to enforce the projection onto the 
S z = subspace. In the following we will use the latter variational wavefunction as 
the starting point for more refined quantum Monte Carlo calculations. In particular 
for J2/J1 = 0.5 we have chosen to work with [3 — 0.4 in Eq. (|75|). The possibility 
to restrict to any total spin projection S z = J^i &i allows one to evaluate the spin 
gap by performing two simulations for S z = and S z = 1. As in the triangular 
case, the potential v(r) used was the same in both subspaces and the variational 
parameter 77 was found by optimizing the energy in the S z = subspace. The latter 
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Table 2. Estimates of the ground-state energy per site and their relative accuracy (in brackets) 
for N = 36 and various values of the J2/J1 ratio. VMC: variational Monte Carlo. VMCLS: 
variational Monte Carlo with LS wavefunction. SRLS: GFMCSR with LS wavefunction and p = 8 
(see text). 



J2/J1 


VMC 


VMCLS 


SRLS 


Exact 


0.0 


-0.6695 (1.4%) 


-0.6756 (0.5%) 


-0.6789 (0.00%) 


-0.67887 


0.1 


-0.6284 (1.5%) 


-0.6349 (0.5%) 


-0.6379 (0.03%) 


-0.63810 


0.2 


-0.5884 (1.8%) 


-0.5952 (0.6%) 


-0.5988 (0.04%) 


-0.59905 


0.3 


-0.5495 (2.3%) 


-0.5574 (0.9%) 


-0.5619 (0.10%) 


-0.56246 


0.4 


-0.5120 (3.3%) 


-0.5237 (1.1%) 


-0.5289 (0.16%) 


-0.52974 


0.5 


-0.4783 (5.1%) 


-0.4916 (2.4%) 


-0.5022 (0.32%) 


-0.50381 



spin-wave variational wavefunction (see Tab. ||) provides a rather good estimate of 
the ground-state energy for J2/J1 < 0.3. Instead such accuracy abruptly decreases 
instead in the regime of strong frustration, suggesting a change in the nature of the 
ground state. 

Within the finite-size spin-wave theory, we can also gain information about the 
low-lying excited states. As shown in Sec. 3.1.2, in order to stabilize a low-energy 
total spin excitation S a magnetic field h in the ^-direction must be added to the 
spin Hamiltonian. Keeping into account that, for small fields, the classical minimum 
energy configuration is the Neel order canted by an angle 6 along the direction of 
h (with sin# = h/2J\z) 1 the finite-size spin-wave expansion is straightforward and 
leads to a linearized Hamiltonian as in Eq. ( fio|) with the following field-dependent 
coefficients 



A^ = l + / 3(«5 k -l)+ 7k (^A-) 2 , B£ 



= -7k 



1 - 



2Jiz 



(76) 



and with a singular part given by 

J\SZ 



n 



SM 



A'o + ^f{S z - Ns sinf 



(77) 



favoring a value of S z (in the original spin representation) consistent with the ap- 
plied field, at the classical level. 

The total spin S = N (S?) of the excitation can be related to the magnetic field 
h by means of the Hellmann-Feynman theorem 




k^0 



(78) 



where 



E(h) = E ct - (shf 



N 



J\SZ 



(79) 
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Figure 12: Spin- wave uniform susceptibility as a function of J2/J1 for s = 1/2. 
The inset displays the low-energy spin-wave spectrum for N — 144 for two values 
of the J2/J1 ratio as a function of |S| 2 = S(S + 1). 



and e£ = y (Afy 2 — {B^) 2 . The final step in order to evaluate the energy spectrum 
E(S) is to perform a Legendre transformation E(S) = E(h) + hsS. 
Finally, the spin-wave uniform susceptibility, 

Xsw = -l/Nd 2 E(h)/dh 2 \ h=0 , 

is, at the leading order in l/.s, 

where Xd = 1/2 J\z. 

4.2. Transition to a quantum disordered state induced by frustration 



4.2.1. Spin-wave susceptibilities and low-energy spectra 

For the unfrustrated Heisenberg model, even for s = 1/2, the spin-wave pre- 
dictions are very accurate as far as the e nergy, the order parameter and the spin 
uniform susceptibility are concerned. t2rt3 Turning on J2 the model is increasingly 
frustrated and one can expect the spin-wave theory to remain accurate only in the 
region where the nature of the order parameter is the same as in the classical case 
(S — » 00). Within this analytical approach we can therefore detect a non-magnetic 
phase by looking for the breakdown of the spin-wave expansion. 

As pointed out by Zhong and Sorella,0 for moderate frustration (J2/J1 < 0.2) 
the linear spin- wave predictions on finite sizes are quite accurate for both the energy 
and the antiferromagnetic order parameter. Moreover, in this regime, the second 
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Figure 13: Low-energy excitation spectra as a function of |S| 2 = S(S + 1) for 
J2/J1 — 0.2, s — 1/2 and N = 16,36,100: spin-wave (full circles and continuous 
line), exact (empty circles and dashed lines), GFMCSR with p = 8 (empty squares 
and dashed lines). 



order correctionEJ leads to an almost exact result. For large J2/J1, instead, the sec- 
ond order term does not improve the first order estimate and a possible breakdown 
of the spin- wave expansion may occur even well below the classical transition point 
J2/J1 = 0.5. In particular at the leading order in 1/s. and for s = 1/2 the order 
parameter vanishes at a critical value (J2/ ' J\) c — 0.38.E3EJ 

Analogously, a breakdown of the spin- wave expansion can be evidenced from the 
vanishing of the uniform susceptibility, which is always finite when there is long- 
range Neel order in the thermodynamic limit and vanishes instead in presence of 
a finite triplet gap (see Sec. 2.2). As it is shown in Fig. 12, the classical uniform 
susceptibility is strongly rcnormalizcd by the quantum fluctuations for s = 1/2 at 
the spin- wave level (jso|). As expected, increasing the frustration such reduction is 
enhanced and leads eventually to the vanishing of the susceptibility for J^j J\ — 0.35. 
The structure of the finite-size spin- wave excitation spectrum below and above this 
critical point is very different (see the inset of Fig. |l2|) with an evident breakdown 
of the quantum top law (^), as well as of the spin- wave approximation scheme, 
in the non-magnetic phase. Below the critical point, instead, the spin-wave theory 
reproduces remarkably well the exact and GFMCSR results for the low-energy part 
of the spectrum in the whole range of sizes (see Fig. |l3| ). Furthermore, as already 
observed for the triangular antiferromagnet, increasing the size, the slope of E(S) 
vs S(S + 1) decreases and gives rise to the collapse of a macroscopic number of 
states with different S on the ground state as N — > oo: i.e., a ground state with a 
broken SU(2) symmetry. 



4.2.2. Size-scaling of the spin gap 



The spin-wave prediction for the occurrence of a non-magnetic region in the 
phase diagram of the J\—J% Heisenberg model is confirmed by our results for the 
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Figure 14: Size scaling of the energy gap to the first S — 1 spin excitation obtained 
with the GFMCSR technique for J 2 /Ji = 0.38 (full triangles), 0.45 (full squares) 
and 0.50 (full circles). Data for the unfrustrated (J2 =0) Heisenberg model taken 
from Ref. El, are also shown for comparison (empty circles). Lines are weighted 
quadratic fits of the data. 



spin triplet gap obtained using the GFMCSR. The latter calculation, which extends 
the recent one by Sorellajlll has been performed using ( [74] ) as guiding wavefunction 
and including in the SR conditions the energy, all Sf independent by symmetry, 
and the antifcrromagnetic order parameter. The latter, though not improving the 
accuracy of the calculation, allows a very stable and reliable simulationfor large p. 
The new results, extended up to N = 144, confirm the previous findingsEH of a finite 
spin gap in the thermodynamic limit for J2/J1 > 0.40 (Fig. [u]). Remarkably, these 
results are not an artifact of the chosen guiding wavefunction: in fact, unlike the FN 
approximation, the GFMCSR is able to detect a finite gap in the thermodynamic 
limit by starting from a spin-wave like wavefunction (|74|) which is Neel ordered 
and therefore gapless. This behavior is very different from the one observed in 
the case of the s = 1/2 Heisenberg antiferromagnet on the triangular lattice (see 
Fig. |Io| ) where, with the same numerical scheme, a similar guiding wavefunction 
and a comparable accuracy we obtained a gapless excitation spectrum. Therefore 
the existence of a gapped phase in the regime of strong frustration is likely to be a 
genuine feature of the J1—J2 Heisenberg model. 

4.3. The nature of the non-magnetic phase 

In principle either a RVB crystal, with some broken spatial symmetry, or a 
homogeneous spin liquid is compatible with a triplet gap in the excitation spectrum. 
Among the dimerized phases proposed in the literature, the so-called columnar and 
plaquette RVB are the states which are the most likely candidates. These kind of 
states can be thought of as a collection of valence bond states | • • ) = |Ti) — lit) 
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Figure 15: Columnar (a) and plaquette (b) RVB states. 



between neighboring sites arranged in the patterns shown in Fig. [15], where the 
plaquettes in (b) are the following rotationally invariant superpositions: 



□ 



+ 



I I 



Both the columnar and the plaquette states break the translation invariance along 
the x and y directions, but only the latter preserves the symmetry of interchange 
of the two axes. 

Read and SachdevEll with a field-theoretic large-iV expansion, were the first 
to conjecture that quantum fluctuations and a next-nearest-neighbor frustrating 
interaction could drive the ground state of the square lattice antiferromagnet into 



have supported over the 
with a study that combines 



a columnar RVB state and series expansion studie; 
years this prediction. Recently, Kotov and co-workerall 
an analytic effective Hamiltonian approach, extended dimer expansions and exact 
diagonalizations have presented a body of evidences that has been interpreted as 
supporting the columnar scenario. Finally, using the GFMCSR with a Density 
Matrix ilenormalization Group guiding wavefunction, du Croo de Jongh and co- 
workers^!! have proposed a ground state with intermediate properties between the 
plaquette and columnar RVB. 



4.3.1. The method of generalized susceptibilities 

In order to better characterize the nature of the ground state in the gapped 
phase, we have checked the occurrence of some kind of crystalline order, by calcu- 
lating the response of the system to o p erators breaking the most important lattice 
symmetries. As suggested in Refs. EjOO, and also shown in Sec. 2.2, the occurence 
of some kind of crystalline order in the thermodynamic ground state can be checked 
by adding to the Hamiltonian (|^) a term —SO, where O is an operator that breaks 
some symmetry of TL. In fact, if true long-range order exists in the thermodynamic 
ground state, the finite-size susceptibility xo = (0}s/NS has to diverge with the 
system size and, in particular, it is bounded from below by the system volume 
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Figure 16: Exact results for the J1 — J2 chain: xt(^) associated to the operator 
Ot (breaking the translational invariance) for J2/J1 = 0.2 (a) and J2/J1 = 0.4 (b). 
Data are shown for TV = 12, 14, 16, 20, 24, 30 for increasing values of xt($)- Lines 
are guides for the eye. 



squared [Eq. (|27|)]. Thus susceptibilities are in principle a very sensitive tool - 
much more than the square of the order parameter - for detecting the occurrence 
of long-range order. 

Within a numerical technique, the susceptibility xo — d 2 e(S)/dS 2 \s=o can be 
calculated with only energy measurements by computing the ground-state energy 
per site in presence of the perturbation for few values of S and by estimating nu- 
merically the limit 

Xo = lim XoW = -«^). (81) 

As we have tested in the one dimensional J\ — J2 model, the numerical study 
of long-range order by means of x(^) is very effective and reliable. Here a quantum 
critical point at J2/J1 — 0.2412 separating a gapless spin-fluid phase from a gapped 
dimerized ground state (which is two-fold degenerate and adiabatically conne cted to 
the Majumdar-Ghosh exact solution for J2/J1 = 0.5) is rather well accepted. EjOlJ 
As shown in Fig. E^, the response of the system to the perturbation SOt, with 

(h ^V ? "S, -S.,.„. . (82) 

i 

breaking the translation invariance with momentum k = ir, is very different below 
and above the dimer-fluid transition point. However it is extremely important 
to perform very accurate calculations at small S to detect the divergence of the 
susceptibilities for large system sizes. 



4.3.2. Stability of Plaquette vs Columnar RVB 
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As also suggested in a recent paper by Singh and co-workers, E3 the appearance 
of a columnar state can be probed by using as order parameter the operator 

Oc = ^ (Si ■ S i+X ~ Si ■ Si +y ) , (83) 

i 

where x = (1, 0), y = (0, 1). As shown in Fig. [I?], the exact diagonalization results 
for N — 16 and N — 36 indicate that the susceptibility associated with this kind 
of symmetry breaking, xci decreases with the system size. In order to exclude 
an anomalous size scaling we have extended the calculation up to iV = 64. Our 
quantum Monte Carlo results, which reproduce quite well the ED data, rule out 
clearly the columnar dimerization. 

The above result is in disagreement with the conclusions of several series expan- 
sion studies 00 However, as stated in Ref. E3, the series for \C are very irregular 
and do not allow a meaningful extrapolation to the exact result. In our calculation 
instead, even the ED results for N < 36, are already conclusive. 

Having established that the columnar susceptibility is bounded, it is now im- 
portant to study the response of the J\ — J2 model to a small field coupled to the 
perturbation 

6 T = J2e iQori Si-S i+x , (84) 

i 

with Qo = (7r, 0), explicitly breaking the translation invariance of the Hamiltonian. 
The evaluation of xt, with a reasonable accuracy, is a much more difficult task. In 
fact in this case the ED values of the susceptibility for N = 16 and N = 32 increase 
with the size and much more effort is then required to distinguish if this behavior 
corresponds to a spontaneous symmetry breaking in the thermodynamic limit. As 
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Figure 18: xt(S) associated to 6 T for J 2 /Ji = 0.5, N = 32 (a), N = 64 (b) and 
N = 100 (c): FN (empty squares), GFMCSR (full squares), FN with LS (empty 
circles), GFMCSR with LS (full circles), exact (empty triangles). 



it is shown in Fig. |l8| (a), the FN technique, starting from a guiding wavefunction 
without dimer order, is not able to reproduce the actual response of the system 
to Ot, even on small sizes. The GFMCSR technique allows us to get an estimate 
of the susceptibility which is a factor of three more accurate, but not satisfactory 
enough. In order to improve on this estimate, we have attempted to include in the 
SR conditions many other, reasonably simple, correlation functions (such as the 
spin-spin correlation functions Sj • Sj for | r*i — | > V2), but without obtaining a 
sizable change of the estimate of xt- In fact, the most effective SR conditions are 
those obtained with operators more directly related to the Hamiltonian.E3o 

After many unsuccessful attempts, we have realized that it is much simpler and 
straightforward to improve the accuracy of the guiding wavefunction itself. This can 
be obtained by applying a generalized Lanczos operator (1 + aTL) to the variational 
wavefunction \4>g)i where a is a variational parameter. This defines the so-called 
one Lanczos step (LS) wavefunction.L3 In the present model by using the LS wave- 
function, a clear improvement (by about a factor of 3) on the variational estimate 
of the ground-state energy is obtained at all strengths of frustration (see Tab. |^). 
With this starting point the GFMCSR provides an estimate of the ground-state en- 
ergy which is basically exact for moderate frustration and remarkable accurate for 
Ji/Ji — 0.4 and 0.5. More importantly, as shown in Fig. 18 (a), the LS wavefunction 
allows a much better estimate of the susceptibility. This calculation was obtained 
by including in the SR conditions the energy, the spin-spin correlation functions up 
to next-nearest- neighbors, distinguishing also S^S* and (SfSJ + SfSj) (p = 4). 
The mixed averages of these correlation functions can be computed over both the 
wavefunction \iPg) and the LS wavefunction (1 + aH)\^G) during the same Monte 
Carlo simulation. Thus with a LS wavefunction one can also easily double the num- 
ber of constraints that are effective to improve the accuracy of the method (p = 8). 
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In this case, we have tested that it is irrelevant to add further long-range correlation 
functions in the SR conditions even for large size. 

By increasing the size [see Figs. Q, the response of the system is very strongly 
enhanced, in very close analogy with the one dimensional model in the dimerized 
phase [see Fig. |l6| ( b)]. This is obtained only with the GFMCSR technique, since 
as shown in Fig. |l8|, the combination of FN and Lanczos step alone, is not capable 
of detecting these strongly enhanced correlations. For N — 100 the GFMCSR 
increases by more than one order of magnitude the response of the system to the 
dimerizing field. Of course, for a definite conclusion, one should check whether the 
susceptibility diverges as the volume squared, as implied by Eq. ( |27|) . However, 
in order to obtain quantitatively reliable zero-field extrapolations (pl|), the limit 
of very small fields has to be reached. This is in general possible within exact 
diagonalization (see Fig. |l6|) but it is rather difficult within a stochastic technique 
like the GFMC which is always affected by a statistical error. 



5. Conclusions 

In this paper we have studied the interplay between frustration and zero-point 
quantum fluctuations in the ground state of the triangular and J\ — Ji Heisen- 
berg antiferromagnets. These frustrated systems are the simplest examples of two- 
dimensional spin models in which quantum effects may be strong enough to destroy 
the classical Neel order, thus stabilizing a ground state with symmetries and cor- 
relations different from their classical counterparts. For this reason, in the last few 
years, they have attracted much theoretical interest even if a general consensus on 
the nature of their ground state has not yet been achieved. With this work, by 



using several tec 
Reconfiguration, 



miqucs including the Green function Monte Carlo with Stochastic 
l3@ a quantum Monte Carlo method recently developed to keep 
under control the sign problem, we have put on firmer grounds the conclusions on 
the ground-state properties of these frustrated models. 

Despite the fact that the spin-half Heisenberg antiferromagnet on the triangular 
lattice was the first historical candidate for a non-magnetic ground statej2-!'E3 all 
our results point toward the existence of zero-temperature long-range Neel order. 
In fact, our quantum Monte Carlo simulations provide robust evidences for a gap- 
less spectrum and for a value of the order parameter that, although reduced (by 
about 59%) with respect to the classical case, remains finite in the thermodynamic 
limit. This_js partially in agreement with the conclusions of finite-temperature 
calculations^ suggesting a ground state with a small but nonzero long-range anti- 
ferromagnetic order and with series expansions studiesa indicating the triangular 
antiferromagnet to be likely ordered but close to a critical point. However, in our 
simulation, which to our knowledge represents the first attempt to perform a sys- 
tematic finite-size scaling analysis, the value of the thermodynamic order parameter 
is sizeable, indicating the presence of stable long-range order. Moreover, the accu- 
racy of the finite-size spin-wave predictions indicates that the spin-wave theory is 
a reliable analytical approximation to describe the ground-state properties of the 
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present model. In particular, the effectiveness of the spin- wave theory in reproduc- 
ing on finite sizes the low-energy excitation spectrum provides further support to 
the existence of long-range Neel order in the ground state, suggesting also that the 
value of the uniform spin susceptibility should be very close to the spin- wave result. 
We believe that our results, together with the clear indications recently provided 
by Bernu and co-workersES with a symmetry analysis of the low-energy excitation 
spectra, finally solve the issue of the ordered nature of the ground state of the 
Heisenberg antiferromagnet on the triangular lattice. 

The effects of quantum fluctuations are more remarkable in the J\—Ji Heisen- 
berg model, where the combined effect of frustration and zero-point motion in- 
terferes with the mechanism of spontaneously broken symmetry, giving rise to a 
non-magnetic ground state of purely quantum-mechanical nature. In fact, our spin- 
wave, exact diagonalization, and quantum Monte Carlo results indicate that quan- 
tum fluctuations are able to melt the antiferromagnetic long-range order in the 
regime of strong frustration, driving the ground state into a quantum disordered 
phase at J2I J\ — 0.4. In addition, with Lanczos and quantum Monte Carlo calcula- 
tions we have studied the susceptibilities for the most important crystal symmetry 
breaking operators. Ourresults, while casting serious doubt on the conclusions of 
series expansion studiesj£MI] indicate the plaquette RVB, with spontaneously bro- 
ken translation symmetry and no broken rotation symmetry, as a more plausible 
ground state in the non-magnetic phase. However further investigation is needed to 
clarify the nature of the ground state in thenon-magnetic phase and a more refined 
numerical investigation is now in progressEl]. In the ordered phase, instead, simi- 
larly to the triangular case, we find a remarkable agreement between the spin-wave 
low-energy excitation spectrum and the exact and quantum Monte Carlo results. 
This suggests that the value of the uniform spin susceptibility should be very close 
to the spin-wave prediction up to J2/J1 — 0.30. 

Our results could be also verified experimentally on the novel realizations of these 
frustrated models, like the triangular K/Si(lll):B interfaced and the Li2VOSi04, 
Li 2 VOGe04 compoundsja quite recently argued to be well described by a spin- 
half J1 — J2 Heisenberg model on the square lattice. Forthcoming measurements 
under pressures could also allow one to tune the J2 / J\ ratio and to investigate the 
properties of these systems in various regimes of frustration. 
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